module suty
  !
  ! Fortran Interface for sutyA.s
  !
  !
  ! LICENSE AGREEMENT:
  !
  ! You agree that we don't agree to anything.
  !
  !
  ! DISCLAIMER OF WARRANTIES:
  !
  ! This software is provided to you "AS IS," and the authors and
  ! contributors disclaim any warranty or liability obligations to you
  ! of any kind, implied or express, including but not limited to the
  ! warranties of merchantability, fitness for a particular purpose,
  ! and non-infringement.
  ! 
  implicit none
  !
  interface i2a
     module procedure i2a_dw
     module procedure i2a_qw
  end interface

  !
  interface a2i
     module procedure a2i_L
     module procedure a2i_I
  end interface

  !
  interface a2f
     module procedure a2f_L
     module procedure a2f_R
  end interface

  !
  private partition_r8
  private partitionP_i
  private quickSortP_i
  private pivot_r8
  private pivot_i
  interface 
     subroutine unpack2(n,pxyz,pout,pw)
       integer,intent(in),value::n,pxyz,pout,pw
     end subroutine unpack2
  end interface
  interface 
     subroutine cp2d(pa,pb,n)
       integer,intent(in),value::pa,pb,n
     end subroutine cp2d
  end interface
  interface 
     subroutine cle2d(pa,n)
       integer,intent(in),value::pa,n
     end subroutine cle2d
  end interface
  interface 
     integer function mcpf(dst,src,len)
       integer,intent(in),value::dst,src,len
     end function mcpf
     subroutine mcps(dst,src,len)
       integer,intent(in),value::dst,src,len
     end subroutine mcps
  end interface
  interface 
     integer function msetf(dst,src,len)
       integer,intent(in),value::dst,src,len
     end function msetf
  end interface
  interface 
     subroutine msets(dst,src,len)
       integer,intent(in),value::dst,src,len
     end subroutine msets
  end interface

  !
  interface kthsmallest
     subroutine kthsmallest(a,n,k)
       integer,intent(in),value::a,n,k
     end subroutine kthsmallest
  end interface kthsmallest

contains

  !
  recursive subroutine quickSort_r8(a,i,j)
    ! sorting of real 
    integer,intent(in)::i,j
    real*8,intent(inout)::a(:)
    integer p,k
    real*8 x
    if(i==j) return
    p=pivot_r8(a,i,j)
    if(p/=-1) then
       x=a(p)
       k=partition_r8(a,i,j,x)
       call quickSort_r8(a,i,k-1)
       call quickSort_r8(a,k,j)
    end if
  end subroutine quickSort_r8
  !
  subroutine persort_i(n,a,p)
    ! sorting of integer with permutation
    integer,intent(in)::n
    integer,intent(inout)::a(n)
    integer,intent(out)::p(n)
    integer i
    do i=1,n
       p(i)=i
    end do
    call quickSortP_i(a,1,n,p)
  end subroutine persort_i
  !  
  recursive subroutine quickSortP_i(a,i,j,per)
    ! this will cause bus error with -O option for large number of
    ! iterations
    integer,intent(inout)::a(:)
    integer,intent(in)::i,j
    integer,intent(inout)::per(:)
    integer p,k,x
    if(i==j) return
    p=pivot_i(a,i,j)
    if(p/=-1) then
       x=a(p)
       k=partitionP_i(a,i,j,x,per)
       call quickSortP_i(a,i,k-1,per)
       call quickSortP_i(a,k,  j,per)
    end if
  end subroutine quickSortP_i
  !
  integer function partition_r8(a,i,j,x)
    real*8,intent(inout)::a(:)
    integer,intent(in)::i,j
    real*8,intent(in)::x
    integer l,r
    real*8 t
    l=i
    r=j
    do while(l<=r)
       do while(l<j.and.a(l)<x)
          l=l+1
       enddo
       do while(r>i.and.a(r)>=x) 
          r=r-1
       enddo
       if(l>r) exit
       t=a(l)
       a(l)=a(r)
       a(r)=t
       l=l+1
       r=r-1
    end do
    partition_r8=l
  end function partition_r8
  !
  integer function partitionP_i(a,i,j,x,per)
    integer,intent(inout)::a(:)
    integer,intent(in)::i,j
    integer,intent(in)::x
    integer,intent(inout)::per(:)
    integer l,r,k,t
    l=i
    r=j
    do while(l<=r)
       do while(l<j.and.a(l)<x)
          l=l+1
       enddo
       do while(r>i.and.a(r)>=x) 
          r=r-1
       enddo
       if(l>r) exit
       t=a(l)
       a(l)=a(r)
       a(r)=t
       k=per(l)
       per(l)=per(r)
       per(r)=k
       l=l+1
       r=r-1
    end do
    partitionP_i=l
  end function partitionP_i
  !
  integer function pivot_r8(a,i,j)
    ! internal use only
    integer,intent(in)::i,j
    real*8,intent(in)::a(:)
    integer k
    k=i+1
    do while(k<=j)
       if(a(i)/=a(k)) exit
       k=k+1
    end do
    if(k>j) then
       pivot_r8=-1
    elseif(a(i)>=a(k)) then
       pivot_r8=i
    else
       pivot_r8=k
    endif
  end function pivot_r8
  !
  integer function pivot_i(a,i,j)
    ! internal use only
    integer,intent(in)::i,j
    integer,intent(in)::a(:)
    integer k
    IF(I+1<J) THEN 
       k=i+1
       do while(k<j)
          if(a(i)/=a(k)) exit
          k=k+1
       end do
    ELSE
       PIVOT_I=-1
       RETURN
    END IF
    if(k>j) then
       pivot_i=-1
    elseif(a(i)>=a(k)) then
       pivot_i=i
    else
       pivot_i=k
    endif
  end function pivot_i
  !
  character*32 function i2a_dw(int,digit,hex,desc)
    integer*4,intent(in)::int
    integer,intent(in),optional::digit
    logical,intent(in),optional::hex
    character*(*),intent(in),optional::desc
    i2a_dw=i2a_qw(int8(int),digit,hex,desc)
  end function i2a_dw
  character*32 function i2a_qw(int,digit,hex,desc)
    integer*8,intent(in)::int
    integer,intent(in),optional::digit
    logical,intent(in),optional::hex
    character*(*),intent(in),optional::desc
    integer rep
    if(int==0.and.present(desc)) then
       i2a_qw=trim(desc)
       return
    end if
    if(present(hex)) then
       write(i2a_qw,"(z16)") int
    else
       write(i2a_qw,*) int
    end if
    i2a_qw=adjustl(i2a_qw)
    if(present(digit)) then
       rep=digit-len_trim(i2a_qw)
       if(rep>0) i2a_qw=repeat("0",rep)//trim(i2a_qw)
    end if
  end function i2a_qw
  !
  !
  logical function a2i_L(a,i)
    character*(*),intent(in)::a
    integer,intent(out)::i
    integer x,istat
    read(a,*,iostat=istat) x
    if(istat==0) then
       i=x
       a2i_L=.true.
    else
       a2i_L=.false.
    end if
  end function a2i_L
  integer function a2i_I(str)
    character*(*),intent(in)::str
    integer i,istat
    a2i_I=0
    read(str,*,iostat=istat) i    
    if(istat==0) a2i_I=i
  end function a2i_I
  !
  logical function a2f_L(str,x)
    character*(*),intent(in)::str
    integer istat
    real f,x
    read(str,*,iostat=istat) f
    if(istat==0) then
       x=f
       a2f_L=.true.
    else
       a2f_L=.false.
    end if
  end function a2f_L
  real function a2f_R(str)
    character*(*),intent(in)::str
    integer istat
    real f
    a2f_R=0
    read(str,*,iostat=istat) f
    if(istat==0) a2f_R=f
  end function a2f_R
  !
  character*32 function f2a(f,fdigit)
    real,intent(in)::f
    integer,intent(in),optional::fdigit
    integer i,j,k,p
    character*32 ex
    write(f2a,*) f
    f2a=adjustl(f2a)
    if(present(fdigit)) then
       ex=""
       j=index(f2a,"E") ! the first exp
       if(j==0) j=index(f2a,"e")
       if(j/=0) then
          ex=trim(f2a(j:))
          p=0
          do i=2,len_trim(ex)
             ! remove leading zeros from exp
             if(ex(i:i)=="0".and.p==0) then
                p=i
             else if(ex(i:i)/="0".and.p/=0) then
                ex(p:)=ex(i:)
             end if
          end do
          k=j-1
          f2a(j:)=""
       end if
       k=len_trim(f2a)  ! the last digit
       i=index(f2a,".") ! dot 
       if(i>0) then
          if(fdigit<0) then
             ! remove trailing zeros
             do i=k-1,i+1,-1
                if(f2a(i:i)/="0") exit
                f2a(i:)=""
             end do
          else
             i=i+fdigit+1
             if(i<=32) f2a(i:)=""          
          end if
       else if(fdigit>0.and.len_trim(f2a)+fdigit+1<=32) then
          f2a=trim(f2a)//"."//repeat("0",fdigit)        
       end if
       i=len_trim(f2a)
       if(f2a(i:i)==".") then
          f2a(i+1:i+1)="0"
       else
          if(index(f2a,".")==0) f2a=trim(f2a)//".0" ! always real
       end if
       f2a=trim(f2a)//trim(ex)
    end if
  end function f2a
  !
  subroutine lower(str)
    character*(*),intent(inout)::str
    integer*1 c
    integer i
    pointer(p,c)
    do i=0,len_trim(str)-1
       p=loc(str)+i
       if(c>=z'41'.and.c<=z'5a') c=int(ior(c,z'20'),1)
    end do
  end subroutine lower
  !
  subroutine strapp(a,b)
    character*(*),intent(inout)::a
    character*(*),intent(in)::b
    a=adjustl(trim(a)//" "//trim(b))
  end subroutine strapp

  !
  real*4 function median(n,a)
    integer n
    real*4,intent(inout)::a(n)
    integer k
    median=0.0_4
    if(n<2) return
    if(iand(n,1)==1) then
       k=(n+1)/2
    else
       k=n/2
    end if
    
    call kthsmallest(loc(a),n,k)
    median=a(k)

  end function median

end module suty
